“Calhoun 


Institutional Archive of the Naval Postgraduate School 





Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


2008-03 


Autonomous underway replenishment at sea 
for Riverine operations 


Addison, William F. 


Monterey, California. Naval Postgraduate School 
http://hdl.handle.net/10945/4278 


This publication is a work of the U.S. Government as defined in Title 17, United 
States Code, Section 101. Copyright protection is not available for this work in the 
United States. 


Downloaded from NPS Archive: Calhoun 


Calhoun is the Naval Postgraduate School's public access digital repository for 


\§ D U DL EY research materials and institutional publications created by the NPS community. 
«iis Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


NNN KNOX appointed -- and published -- scholarly author. 


LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 


http://www.nps.edu/library Monterey, California USA 93943 





NAVAL 
POSTGRADUATE 
SCHOOL 


MONTEREY, CALIFORNIA 


THESIS 


AUTONOMOUS UNDERWAY REPLENISHMENT AT SEA FOR 
RIVERINE OPERATIONS 


by 





William F. Addison 











March 2008 


Thesis Advisors: Fotis Papoulias 
Oleg Yakimenko 





Approved for public release; distribution is unlimited 





TH 


IS PAGE 








NTENT 





ONALLY LEFT BLANK 





REPORT DOCUMENTATION PAGE Form Approved OMB No. 0704-0188 


Public reporting burden for this collection of information is estimated to average 1 hour per 
response, including the time for reviewing instruction, searching existing data sources, gathering 
and maintaining the data needed, and completing and reviewing the collection of information. Send 
comments regarding this burden estimate or any other aspect of this collection of information, 
including suggestions for reducing this burden, to Washington headquarters Services, Directorate 
for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 
22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188) 
Washington DC 20503. 


1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 
March 2008 Master’s Thesis 


4. TITLE AND SUBTITLE Autonomous Underway 5. FUNDING NUMBERS 
Replenishment at Sea for Riverine Operations 
6. AUTHOR(S) William F. Addison II 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS (ES) 8. PERFORMING ORGANIZATION 
Naval Postgraduate School REPORT NUMBER 

onterey, CA 93943-5000 

9. SPONSORING /MONITORING AGENCY NAME(S) AND 10. SPONSORING/MONITORING 
ADDRESS (ES) AGENCY REPORT NUMBER 
N/A 











11. SUPPLEMENTARY NOTES The views expressed in this thesis are those of the author and 
do not reflect the official policy or position of the Department of Defense or the U.S. 
Government. 


12a. DISTRIBUTION / AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE 
Approved for public release; distribution is unlimited 


13. ABSTRACT (maximum 200 words) 

Currently the United States Navy is making a small footprint in the world’s littoral 
regions with the help of the United States Marine Corps. In Iraq, the Marine Corps is 
actively conducting Riverine operations, however they are overly tasked and in need of 
permanent replacement by the United States Navy. In order to alleviate the Marine 
Corps, the Naval Expeditionary Combat Command with its Riverine Squadrons will soon 
tak over thes Riverin operational commitments in order to reestablish supremacy 
throughout the Riverine environment. With this in mind, the Chief of Naval 
Operations, Center for Naval Analyses requirements, System Engineering Analysis (SEA- 
11) class of 2007 developed a concept of operations (CONOPS) which the Total Ships 
System Engineering (TSSE) class of 2007 used to develop a prototype platform, which met 
all initial design requirements. In order to take full advantage of this prototype 
platform, every effort was taken in order to minimize the number of crew members on 
station at any given time. The purpose of this thesis is to demonstrate the use of th 

direct method, which will allow the Specialized Command and Control Craft (SCCC) to 
conduct a fully autonomous Underway Replenishment at Sea (UNREP) with a standard supply 
vessel. The direct method approach allows for a smooth path is created instead of 
using waypoint navigation. Additionally, this method allows for real-time updates at 
(1Hz). 



































14. SUBJECT TERMS AUV, UUV, robotics, trajectory planning, path 15. NUMBER OF 
planning, rendezvous, real-time, direct methods PAGES 





TA 


16. PRICE CODE 


17. SECURITY 18. SECURITY 19. SECURITY 20. LIMITATION OF 
CLASSIFICATION OF CLASSIFICATION OF THIS CLASSIFICATION OF ABSTRACT 
REPORT PAGE ABSTRACT 
Unclassified Unclassified Unclassified UU 
NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89) 
Prescribed by ANSI Std. 239-18 








TH 


IS PAGE 








NTENT 





ONALLY LEFT BLANK 





nl 


Approved for public release; distribution is unlimited 


AUTONOMOUS UNDERWAY REPLENISHMENT AT SEA FOR RIVERINE 


OPERATIONS 





William F. Addison 
Lieutenant, United States Navy 
B.S., United States Naval Academy, 2001 











Submitted in partial fulfillment of the 
requirements for the degree of 





MASTER OF SCIENCE IN MECHANICAL ENGINEERING 


from the 


NAVAL POSTGRADUATE SCHOOL 














March 2008 
Author: William F. Addison 
Approved by: Fotis Papoulias 


Thesis Advisor 


Oleg Yakimenko 
Thesis Advisor 


Anthony J. Healey 
Chairman, Department of Mechanical 
Astronautical Engineering (MAE) 











Tal 


and 





TH 


IS PAGE 








NTENT 





ONALLY LEFT BLANK 





iv 


ABSTRACT 


Currently the United States Navy is making a small 





footprint in the world’s littoral regions with the help of 














the United States Marine Corps. n Irag, the Marine Corps 





is actively conducting Riverine operations, however they ar 


overly tasked and in need of permanent replacement by the 








United States Navy. In order to alleviate the Marine Corps, 





the Naval Expeditionary Combat Command with its Riverine 
Squadrons will soon take over these Riverine operational 


commitments in order to reestablish supremacy throughout the 








Riverine environment. With this in mind, the Chief of Naval 











Operations, Center for Naval Analyses requirements, System 
Engineering Analysis (SEA-11) class of 2007 developed a 
concept of operations (CONOPS) which the Total Ships System 
Engineering (TSSE) class of 2007 used to develop a prototype 





platform, which met all initial design requirements. In 








order to take full advantage of this prototype platform, 
every effort was taken in order to minimize the number of 
crew members on station at any given time. The purpose of 


this thesis is to demonstrate the use of the direct method, 








which will allow the Specialized Command and Control Craft 


(SCCC) to conduct a fully autonomous Underway Replenishment 





at Sea (UNREP) with a standard supply vessel. The direct 


method approach allows for a smooth path is created instead 





of using waypoint navigation. Additionally, this method 


allows for real-time updates at (1Hz). 
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I. INTRODUCTION 


A. BACKGROUND 


When the end of the “Cold War” with the Soviet Union 





came about, there was a major shift in the United States 


Naval Doctrine. With the United States Navy’s major 





opponent on the high seas eliminated, so to was the threat 





of fighting a major naval battle on the high seas as was the 
threat in years past. After numerous studies and analysis 


of current naval operations and assets, in August 2005 the 





Chief of Naval Operations (CNO) announced that the United 





States Navy would reconstitute a Riverine capability, 
allowing the United States Navy to transition from a blue 
water navy to a force which would be capable of sustaining 
operations in the littoral regions of the world. The CNO’s 


vision called for the resurgence of the brown water Riverine 








Force which is called out in the CNO’s Concept of Operations 


(CONOPS) for the 21°* Century Riverine Force. This document 





calls for the formation of a Naval Expeditionary Combat 








Command, which requires a Riverine force as one of its 
elements. The primary mission for this force is to conduct 
Phase 0 (shaping and stability) operations, to provide 
maritime security and to carry out additional tasks 


specifically related to the Global War on Terrorism. [1] 








Figure l. Blue Water to Brown Water Navy 


Currently the United States Navy is making a small 


footprint in the world’s littoral regions with the help of 











the United States Marine Corps. n Iraq, the Marine Corps 








is actively conducting Riverine operations, however they are 


overly tasked and in need of permanent replacement by the 





United States Navy. In order to alleviate the Marine Corps, 


the Naval Expeditionary Combat Command with its Riverine 





Squadrons will soon take over these Riverine operational 
commitments in order to reestablish supremacy throughout the 


Riverine environment. 


Based on the Chief of Naval Operations, Center for 


Naval Analyses requirements, System Engineering Analysis 





(SEA-11) class of 2007 developed a concept of operations 
(CONOPS) which the Total Ships System Engineering (TSSE) 
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class of 2007 used to develop a prototype platform, which 


met or exceeded all initial design requirements. The first 








step in the development of this new platform, was to conduct 





a Capability Gap Analysis of existing Naval assets both US 





and foreign. It was quickly determined from this analysis, 


that no current ships were capable of fulfilling all of the 








initial requirements requested and a new platform would be 


needed to accomplish the vision of the CNO. Based off of 





these results, a functional Element Decomposition of the 


system requirements was developed and a preliminary design 





was identified. The ultimate design evolved into a multi- 


hulled Specialized Command and Control Craft (SCCC), which 





would utilize three multi-mission craft (MMC) to accomplish 


all mission requirements. 


The existing procedure for conducting Riverine 





operations is to first establish a land forward operating 


base and to then deploy Riverine assets from this land based 





support center to carryout various missions. While Irag has 





shown a land basing system to b ffective, in the future it 


maybe more likely that the Navy will require a sea based 





support structure in order to accomplish its Riverine 








mission. To accomplish this future scenario, the Navy could 
use existing assets; however, this approach limits the 


Navy’s future Riverine footprint due to the limited access 





current assets have in the majority of the rivers of the 
world due to these assets’ slow speeds and deep drafts. To 
structure the United States Riverine forces in such a way as 


to create maximum operational flexibility in the majority of 








the world’s littoral regions, this new platform must be 


produced. 


B. MOTIVATION 


1. Manning 


The Riverine forces will comprise of 800 personnel 
divided among three squadrons. Minus the command structure 
for the Riverine force, each squadron allotted 225 
personnel. Each squadron will consist of three SCCCs and 
nine MMCs. The TSSE Manning Study resulted in the need for 
216 personnel to fully man these ships, with the remaining 
nine comprising the command structure of the squadron. The 
command structure will remain afloat on the Global Fleet 
Station (GFS) in order to manage the logistics of the SCCC. 


The GFS will be removed from the Riverine Area of Operation 





and in Blue water. 





Figure 2. Bridge Layout 


Every effort has been made in the design of the 


Tiberinus Class to minimize the number of crew members on 
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station at any given time. With this in mind, the bridge 
will be the only actively manned space on the ship. As with 
existing naval ships, the bridge has been organized to allow 


the commanding officer to oversee all aspects of navigation; 





however, the bridge on the Tiberinus Class has also 








incorporated the Combat Information Center and Damage 





Control Central, in an effort to centrally locate all 








controlling stations. All bridge watch stations will be 





equipped with touch screen panels which will enable any 
watch stander to reassign their watch station to receive 
additional information from any other watch station and take 
control from their console. Figure 2 shows a purposed 


bridge layout. With the changes from traditional naval 





watch team structure, it will require a crew of 72 personnel 
to fully man the Tiberinus SCCC and the three associated 
MMCs. 


2. General 





Due to the reduced size and complexity of the Tiberinus 
Class, Significant advances in automation of processes and 
procedures had to be achieved in order to allow its reduced 
crew to fully operate the ship in efforts to achieve all 
mission objectives. As a result, all efforts have been made 


to allow the ship’s crew to operate the ship with minimal 





personnel on station. 


Currently there are no options for ships to conduct 





Underway Replenishment at Sea (UNREPS) operations 
autonomously. The ability to accomplish this would allow 
for increased force flexibility and operation. 


Additionally, if not used for fully autonomously UNREPS, 





this technology could be used as visual cueing for complex 
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formation maneuvers, UNREPS, plane guard operations, or even 
pier dockings. This technology could be incorporated with a 


heads up display, which would use standard maneuvers to 








build a database of near-optimal trajectories calculated 


beforehand. These near-optimal trajectories would allow the 




















Officer of the Deck (OOD) to not just mentally visualize the 














command, but this technology would allow the OOD to actually 





see a Simulation of where he or she will end up, thus adding 


to the overall situational awareness. 





Cc. SCOPE 


The scope of this thesis will consist of analytically 
developing a path-planning process which will generate 


trajectories for an UNREP between a standard USNS oiler and 





the SCCC. The first step in achieving this objective will 
be to develop a hydrodynamic model of the SCCC in order to 








utilize the equations of motion in which it will be used to 





Simulate the vessel movement. 





The next step will be to formulat th rendezvous 





trajectories based off of mathematical basis. With this 





information, the factors which will effect the trajectories 


shape can be explained and constraints can be formulated to 





take into account both permissible trajectories and the 





vessel constraints. From here, this information will be used 





to generate the performance index of the vessel. 





After the trajectories are computed, the inverse 


dynamics will then be used to calculate the required states 








of the vessel at each point upon the trajectory path. In 








order to minimize any violations of the optimal parameters, 





the values generated by the trajectory algorithm and the 


inverse dynamics will be used in the performance index. 





D. PROBLEM FORMULATION 
The problem can be summed up as follows: The supply 
vessel will provide a rendezvous point. From this point, as 


an example they will suggest a course of 090 degrees, speed 


13 knots; however, as is often the case, they may need to 





maneuver in order to avoid a contact or to creat a better 





UNREP situation. As the supply vessel maneuvers, they will 





send updates to the SCCC so that it may mak real-tim 








updates in order to achieve station at a lateral separation 





distance of 140 ft, while maintaining C090/S13kts. 


Once this is achieved, a laser range find system will 
be employed to maintain the lateral separation. This is 


illustrated in the figure below. 


Ss i ie ee 
os 7 » “a =— = D> USS SUPPLY 
Fa es 


-- neta =o. 


SCCC 






Updates based 
on course 


Figure 3. Problem Outline 


The figure above shows a pictorial representation of 





the problem as described above. The proposed sequence of 
events is to have the supply vessel communicate with the 
SCCC in order to command the SCCC to proceed to the 
rendezvous point. The SCCC will then compute the necessary 


trajectory to complete the mission and reply with an 





acknowledgement or the SCCC will decline the command 
request. A denial from the SCCC would constitute a 
violation in one of the system constraints and a request by 
the SCCC would then be sent in order to allow the SCCC to 


reach the rendezvous point by either altering the time of 
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+ 
140 ft 
i 





arrival or by requesting a different rendezvous point all 
together. This update would be achieved at a rate of (1 Hz) 


or real-time. 


E. THESIS STRUCTURE 


The intent of this research is to develop a direct 
method control for the SCCC in order to allow for an 


autonomous UNREP. This will be conducted by first 





determining the equations of motion, followed by the 





development and validation of the trajectories, and finally 


by the vessel simulations. 














Chapter will focus on the Tiberinus Class. It will 














focus on the equations of motion for the SCCC and the vessel 























simulation development. Chapter explains the theory and 








equations used for the direct method for rapid prototyping. 











Chapter IV will focus on the development and validation of 


the trajectories for the SCCC. Chapter V will present a 





Simulation for the UNREP between the SCCC and a_ supply 
vessel. Additionally, Chapter V will provide the thesis 


conclusions. 
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II. THE TIBERINUS CLASS 


A. SPECIALIZED COMMAND AND CONTROL CRAFT DESCRIPTION 





The Tiberinus Class Ship has been designed to provide 


command, control and support for Shaping and Stability 





operations. In support of the previous mentioned 
operations, the SCCC will provide Maritime Security and 
carry out additional tasks specifically related to the 


Global War on Terrorism (GWOT) throughout the littoral 








regions of the world. By design, the Tiberinus Class will 





provide a sea-based maritime capability which will enable 


U.S. forces to have an enhanced presence in their areas of 








operation, will maintaining the legitimacy and sovereignty 


of the United States’ ally and coalition partners lands. 





Figure 4. SCCC and MMC 
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This class has been designed to allow for a Riverine 
force to sustain a forward presence within a Riverine 


environment anywhere in the world for an indefinite period 





of time, while maintaining the capability of conducting 
interdiction operations, low intensity combat operations, 


Visit Boarding Search and Seizure (VBSS), maritime security 





operations, and waterborne checkpoints. 


Each Tiberinus Class Ship has been fully designed to be 


independent of each other and will provide all of its own 








hotel services for its embarked personnel. Although the 
Tiberinus Class has primarily be designed in a supporting 


role, each ship has been designed to allow it to carryout 





the same missions as the MMCs in reference to a combat role. 
The characteristics of this class are as shown in the table 


below. 


Characteristics 


(LOA) 
68 ft 
40+ kts 
Ga Se 











(Design) LeSUO sam 
(Maximum) 3,750 nm 


Displacement 
(Design) 








550 LT 





1 H-60 (landing, not 
Aircraft housed) 


Mission Craft 3 JMECs / MMCs 








Table 1. Ship Characteristics 
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B. EQUATIONS OF MOTION 





In designing the model for this vessel, the first step 





is to determine th quations of motion. Since the SCCC and 
the supply vessel will be operating in two dimensions, only 
equations in the horizontal plane will be considered, 


however, this process will be described for a three 





dimensional application for future iterations. This is a 
somewhat common equation development and will be briefly 


described below. 


When dealing with relative motion certain assumptions 





must be established in regards to the motion boundaries. 


For the purposes of this research, it will be assumed that 





the vessel will act as a rigid body, which will enable for 





the calculation of forces and moments on the vessel. 





Additionally, it will be assumed that the Earth’s rotation 








is negligible in regards to acceleration components of the 





vehicle’s center of mass. This will allow for the illusion 





that the vessel will be moving over a stationary plane. 


Finally, it will be assumed that the forces acting on the 








vessel will have their origins in an inertial and/or a 





gravitational prospective. Additionally, the other primary 


forces action on the vessel will be hydrostatic, propulsion, 








thruster, and hydrodynamic forces from lift and drag. 


For vehicles and vessels described in terms of three 





dimensional components, the velocity of these 





vehicles/vessels will be accounted for using six terms given 





as surge, sway, heave, roll, pitch, and finally yaw. For 


the purpose of this research, surge (u) is the vessels 





forward speed; sway (v) is the side slip velocity, heave (w) 
corresponds to a velocity component in the local 2 
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direction, however, its global velocity components do depend 








on the vessels heading, pitch, and roll. These three 


components makeup the body fixed coordinate as shown below. 


R, =T(a@)[wi+vj+wk] (2-1) 





The other three terms are represented of the angular rate of 





rotation of the bodies fixed frame @o as shown below. 
o=[pi+q t+rk] (2-2) 


The vector quantity of these additional three terms, have 








their defined meanings in terms of the vessel motion. The 
vessel “roll rate” is described by the p component in 
equation (2-2). The vessel “pitch rate” is described by the 
gq component in equation (2-2). Finally, the “yaw rate” of 
the vessel is described by the r component as seen in 
equation (2-2). These particular components would be sensed 
by the onboard gyro of the SCCC. All of these components 


combined, account for the overall velocity of the vessel, 





which can be displayed as shown below or as later displayed 














in Chapter 





x=[uvwpqrl (2-3) 


All of the applied external loads for the body 


coordinate components ar represented by the vector 








components of forces which are applied on the body of the 





vessel and the moments which are also applied at the center 


of the body fixed frame as shown below. 


>< 


app 


me 


(2-4) 


dv. 
mS sap, +maxar a +@xv)—f, =| Y,,, 


N 


app 
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I@+ox(I,o)+mip,xv+9,xoxv}—-m, =| M,,, (2-5) 


F(t)=[X ap) Yop (0) Zapp (t) Kapp) Map (1) NayOV (2-6) 





Due to the SCCC being in its first iteration of design, 





the hydrodynamic coefficients were computer simulated in 





order to determine the necessary hydrodynamic coefficients 


for the equations of motion for the vessel. For the purpose 





of this thesis, the constants and hydrodynamic coefficients 
for the SCCC are not included due to the non-trivial nature 


of the task and the unconventional hull shape and type. 





It is often difficult to assess the vessels mass 
moments of inertia about its Center of Gravity (CG), due to 
the CG change with loading and unloading of the vessel. 
These loading and unloading changes are usually symmetric, 


which is ideal when attempting to maintain the vessels 





proper trim and heel under normal static conditions. 
Typically, the mass and angular motion of the vessel is 
described through the mass moment of inertia matrix, which 


is shown below. 
L=|l, f ft (2-7) 


The elements in the above matrix can be determined by the 





following equations below. 


N 
1 = Yd, (9 He): 228) 
i=l 
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y dm (x*+z7) (2-9) 

' Yam (x +y) (2-10) 
bg = Ly = Sd, (xy) (2-11) 
a =) dm,(xz) (2-12) 


N 
1, =1, =->,dm (yz) (2-13) 


i=l 


With the angular velocity vector express as a column vector 


as is shown below, 


sv 


@=|q\ (2-14) 


the angular momentum of the vessel can then b xpressed as: 





H,=I,w (2-15), 


which will allow one to utilize Newton’s second law in order 








to achieve the following equation below. 


2 


dH, @R 
y= Bea pya{m ee (2-15) 





The right hand term in equation (2-15) is representative of 





the moment of the inertial reaction of the sum of the 


external forces acting on the vessel. 


The vertical plane for the equations of motion will 





begin with setting all horizontal plane variables to zero, 
in order to allow only w, q, 96, and z to be the only 
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variables of concern. This, with the addition of a constant 
speed and utilizing small angular changes will allow one to 


utilize a linear system approach. The reduced equations are 





shown below. 





mw, =mU,q + (W —-B)cos6+ 0Z,(t) (2=16) 
1,4 =(Z,B-Z,W) sin + 6Z,(t) (2-17) 
O=4q (2-18) 


One can then manipulate the above equations to create a more 


useable format, which is shown below in the next set of 








equations that will allow the user to conduct matrix 








operations and finally display the state and control 


matrices. 


M x, (t) =a,x,(t)+b,6,(t) (2-19) 


m=Z. =f. <Q 
m,=| —-M, I1,,-M, 0 (2-20) 
0 0 1 
Zs mU,+Z, 0 
a, =|M,, M, 3B =Z.W (2-21) 
0 1 0 


X,(t) = A,x,(t)+ B.O,(t) (2=23) 


xv(t)=[w q Oy (2=24) 
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Av=|-M, 1,-M, 0] |M, M, 2,B-zW| (2-25) 
0 0 1] [0 1 0 
m-Z, -Z, 0] [Z, 
B,=|-M, 1,-M, 0| |M, (2-26) 
0 0 1] [0 


Disregarding the motions in the vertical plane will 








result in the horizontal equations of motions, which are 


shown below. 


mv, =—mr+Y,v.t+y,v,+¥.r+Yr+Y,6,(0) (2227) 
[7 =N,v,+N,v,+N,7+N,r+N,0,() (2-28) 
y=r (2-29) 


As with the previous section in regards to the vertical 


plane, the horizontal dynamic equation is 
m,X,(t) = a,x, (t) +b,0, (0) (2-30) 
where 
x,@)=r yl (2-31) 


Again, as previously shown, this equation converted into a 





more functional form and the state and control matrices 


result as shown below. 
X,(t() =A,x,()+B,6.(t) (2-32) 


-l 


m-Y, -¥Y, O|'[¥, Y¥,-mU, 0 
A,=|-N, .-N, 0| |N, oN, O| (233) 
0 0 1] jo 1 0 


18 


19 





TH 


IS PAGE 








NTENT 





ONALLY LEFT BLANK 





20 


III.DIRECT METHOD FOR RAPID PROTOYPING 


A. HISTORICAL BACKGROUND 


dy General Description of Direct and Indirect Method 





From many years of research, it has been determined 


that the most precise approach in solving an optimal control 





problem is by the variational approach, which is based on 
the Pontryagin’s minimum principle. This indirect approach 


requires the solution of the necessary conditions of 





optimality associated with the infinite dimensional optimal 








control problem rather than optimizing the cost of a finite 





dimensional discretization of the original problem directly. 








Using this method does require advance analytical skill in 














which one must generate numerical solutions of the resulting 


two-point boundary-value problem. The minimum principle is 





used to eliminate the controls in this indirect method 








approach, resulting in a generally nonlinear function of the 


state and co-state variables. This indirect approach allows 





for the generation of benchmark solutions, which will 





generally converge if only excellent initial guesses for the 





non-intuitive con-states are achieved. Additionally, this 


requires the switching structure to be guessed correctly in 





advance. 


A direct method approach allows for rapid trajectory 


prototyping ability. This method uses finite dimensional 





discretization of the optimal control problem to a nonlinear 





programming problem. While the direct method approach does 





not allow for extremely great precision and resolution as 
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that of the indirect method, it does allow for a more 


practical application as its convergence robustness is far 





superior. 


2. History 


The ideal of the direct method approach was first 
developed by Euler in the early 1900’s, when he approached 
the solution of functions as finite sets of variables. This 


approach was achieved by representing acceptable functions 





in the form of infinite power series 
HG Dee a) 
k=0 


or by Fourier series 


“0 S* (a, coskx + b, sin kx), (3-2) 
k k 


y(x) = 7 


or by any series in the form of 


y(x) = » epee (333) 


k=1 


where @,(x) is a given function. Thus, instead of an 


infinite series, the user is only considering a finite 





series, whose solution is simply the function of a set of 





unknown coefficients. 





Ritz too developed a direct method, in which his method 











requires a field problem to be arranged, in which it will be 
used as an integral minimization. Thus, allowing it to be 
used for problems which have variational principles. 


Galerkin obtained approximate solutions to boundary-value 
problems in a simpler way, which is why it is more of a 


universal process. When Galerkin’s method is combined with 
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the interpolation equations of the method of finite 
elements, it allows one to solve both initial and boundary- 


value problems, which Galerkin utilized to solve parabolic 








and elliptic partial differential equations. 








In the 1950’s, aerospac ngineers began to utilize the 


finite element method. Due to vast improvements in 





computing in the following year, this method became 





popularized for numerous numerical simulations of physical 





problems dealing with stress analysis, structural and solid 


mechanics, heat transfer, and fluid mechanics among others. 








Resulting conclusions found that the previous fore mentioned 
methods when applied, will yield approximation for the 
minima from above or the maxima from below, thus, enabling 


the user to utilize them for rapid prototyping of optimal 











solutions or near-optimal solutions. This allows for the 
ability to preset xtrem trajectories and/or controls, 
while allowing for a calculational advantage, while 











providing a near-optimal solution with any varying degree of 





accuracy. 





In the 1960’s, Taranenko applied a similar method to 





that Ritz-Galerkin to flight dynamics involving constraints 





on states and controls. Continuing in his predecessors’ 








methods, he attempted to use continuous, unequivocal and 


differentiable functions which automatically satisfied 





boundary conditions of the function 


Xp KX. 
xX, =X, + —* (r-1,) + ®, (2), i=1,..4 (3-4) 
Tp —To 





as the reference function for the Cartesian coordinates and 





speed. In equation (3-4), tis an argument, while 9,(7) is a 
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continuous, unequivocal, and differentiable function which 


satisfies the boundary conditions®,(7,)=9,(7,)=0. Taranenko 





further suggested the uses of the following equations or any 


of their linear combinations. 


€=7, 





@i(c) =>) A, sinkz (3-5) 


= T,—T 
Di(r)=) A (t-%)*(e-1,)* (3-6) 
k=l 


i (r)=(t-t))™(t—-1,)™ (3-7) 


While there ar no actual limitations, one could use 








any convenient function for their particular task. State 





parameters and controls can then be resolved from the result 


of their inverse flight dynamics. 








Figure 5. Splitting Original Interval 








In order to provide more flexibility without increasing 
the order n(m,m,) , Taranenko recommended the splitting of 
the interval [t>3T,] into pieces, in order to use lower order 


polynomials in order to describe the behavior of the state 
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variables x (i=1,2,3,4). The higher order n(orm,andm,), the 








higher the number of pieces required in the piecewise case, 


thus the closer a near-optimal solution will be to the 





optimal solution. 





Taranenko continued his research with the hopes of 


building a database of trajectories for numerous aircraft, 








in an effort to aid pilots in maneuvering their aircrafts, 








by suggesting the maneuvers optimal trajectory. Due to the 
lack of computing speed of the day, it was discovered that 
the numerous optimization parameters would not allow for 


onboard computation of these trajectories. 


As with most discoveries, this approach was overlooked 


until technology could appropriately catch-up with the 





computing requirements. In 1997, Taranenko’s dream was 














realized by Yakimenko who tested these methods onboard a 








flying laboratory. Yakimenko developed an algorithm which 
computes near-optimal trajectories. These trajectories were 
found to be capable of satisfying all boundary conditions 


prior to the establishment of the trajectories, which 


allowed for decreased calculational time and the ability to 





conduct real time onboard trajectory computation. 


B. MATHEMATICAL DEVELOPMENT 





In order to develop the mathematics for this problem, 


the first step is to formulate an optimal control problem 











which will move the vessel from a starting point (initial 
point) to some final point. If one chooses or is given both 
an initial and final position, which will also serve as the 








problems boundary conditions, a first order polynomial 
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representation of a trajectory linking both positions can be 


developed by using the following formulas: 
X(T) = P, (7) =a, +a,t (3-8) 
y(T)= P(r) =) + br (3 2)) 


As previously mentioned, T is given as any argument. 
One can then solve for any unknown coefficient by 
substituting any value for 7 and solving for these unknown 


coefficients with the following matrix equations: 


1 Ojla x 
|. er 
z 


























Ll 7 Jia, Mee | 
1 OD 
Ole Yo (3-11) 
L1 7 | 5, yy 
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Final 











Initial 





























0 5 1 1 2 
X-axis 


Figure 6. Basic First Order Polynomial 


The straight line represented in the figure above, 





represents the satisfying Of the initial boundary 


conditions. 


While outlining the problem of an UNREP between the 


SCCC and a supply vessel, a straight line trajectory may not 





be the optimal trajectory to accomplish this objective, due 
to changes in course/speed, contact avoidance, and changes 
in rendezvous time. Based on these foreseen events, a 


curved path by be more appropriate in achieving this goal. 





As outlined by other research endeavors and as previously 
mentioned in this thesis, a curved path can be achieved by 
breaking the problem into smaller pieces or as is standard 
naval practice, put in additional waypoints to the final 
position or objective. While this seems simple enough, the 
added waypoints must then account for added boundary 


conditions of that particular leg, and an added time to the 
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next leg must also be taken into account. While this is 





fine for typical naval applications, it does add to the 


overall computational time to achieve the optimal trajectory 
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Figure 7. Curved Trajectory Using Waypoints 





In Figure 7, two waypoints were added to achieve a 


curved trajectory between the initial and final positions. 








It is of interest to note that in order to approximate this 





new path, it must b represented as a polynomial. 
Additionally, an added waypoint is a direct corresponding 


result to an increase in the order of the initial 





polynomial. In this case the corresponding result would be 





a third order polynomial due to the addition of two 





waypoints. If this new curved path is represented in the 


form of a third order polynomial, it will produce a smooth 





curve which is very similar to the outline of the waypoints, 


28 


25 





while maintaining the existing boundary conditions. Using a 
direct method approach, the vessels controls would be taken 
from the produced trajectory and the required states would 


then be produce by using inverse dynamics. 


Assuming the initial boundary conditions remain the 





same, this third order polynomial representation of the 


intended trajectory can be displayed as previously shown 





with only minor changes as shown below: 


3 

HOSP.) => ar (3-13) 
i=0 
3 a 

y(t) =P,(c) = >be (3-14) 
i=0 


10 0 0 |fa] [x10 0 OTe] Tr, 


2 3 2: 3 

1 TZ UT 7, a, x 1 7, tT 7 WD, J 
7 = (3-05) 

2 3 2. 3 
| a ane a | X% ||1 2, ty 7; |B, yy 

2 3 2 3 

x 

lt, t3 73 |L% PA Mee ae Ms Yi 


By taking into account the first derivative of the 





first order polynomial equations, added path curvature 


flexibility will be increased allowing the user to utilize 





boundary conditions to gain the equations for the velocity. 





An example of the process for a single high-order 





polynomial approximation for the coefficients of (A) in a 





“2+2” case are shown below. 
X(T)=X(OV=Xq X(T, =X, (3-16) 
X';(T)) = x'(O) = x'ig X'(T,)=X'y (3=10) 


x" (T)) =x" (0) =x" x"(t,)=x", (3-18) 
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5 
x' (7) =a, +2a,,7 +3a,,0° + 4a,.0° +5a,.0° = hae 


= 2 3 4 5 
X(T) =A +A ,T HG, .T +,0° +440 +4,.7 


5 
= 24" 
= ae 


k=0 


(3-19) 


1 


(3-20) 


k=0 


5 
x" (t) = 2a, + 6a,,0 +12a4,477 + 20a,5r? =) k(k-Nayt*? (3-21) 


C. 





establ 


required to be stationary points. 


Toolbox, 


lish 














k=0 
10 0 0O 0 Ayo 
0 1 0 0 0 a, 
002 0 0 i, 
1 tf, UR ae a 7 a, 
0 1 22, BG; Ar’ St, ay 
0 0 2 6r, 1222 202? |Las | 





HYPOTHETICAL TRAJECTORIES 


In order to create a random trajectory, 


an initial and final po 


one can determine th 








(3-22) 








the user must 


Sition, which will be 


Using Matlab’s Symbolic 


coefficients of the example 


from the previous section as shown below. 





3x,,-9x, 24x, -36x,,  60(%) — 0) 


aio Xin | | Gs 

Gy |=) Xa | | Ga |= 0 18x; 12x, 84x; + 96%, 

Ajo Xj2 a5 0 0 10(x, —X,) 
In order to conduct a single 
approximation using a visual check, 


below in 





the following table. 
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a 
? 0 |r 
—180(%; —%g) 0 t| @-23) 

60% +39) 120%) |] 774 

5 

Re al 
high-order polynomial 


the inputs are set as 


























where, 


The resulting 
meters. 











Table 2. Example Inputs 


£, =I, 2y.ue410 
Xi = {-0.4;-0.1;0.2;0.5} 


plots are the visual confirmation check 
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25 
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Figure 8. Visual Check 
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in 





a 


above, 


mentioned techniques 


before 





Utilizing th 


random trajectory was established using the SCCC concepts of 


operation to establish the conditions for this trajectory. 














This random trajectory is shown below. 


X-axis (x1) (m) 


Random Trajectory Generation 


Figure 9. 
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IV. TRAJECTORY DEVELOPMENT AND VALIDATION 


A. INDENTIFIYING CONTROL PARAMETERS & CONSTRAINTS 


The SCCC UNREP with a supply vessel is a problem 





centered on navigation with the primary controlled 
parameters of concern being course and _= speed, which 
essentially makes this a two dimensional problem. For the 
propose of this thesis, the two vessels in question will be 
in open ocean not constrained by draft, however, the body of 
water may have above water obstructions such as oil 


platforms. 


Viewing this problem, the constraints which come to 


mind are the rudder/speed rule of 15/15 for a total of 30 





overall. Meaning, one could use a speed of 20 knots and a 
rudder of 10 knots for a total of 30. An additional 
constraint will deal with limiting speed and/or rudder 








commands when within approximately 350 feet of the supply 


vessel as to limit the possibility of collision. 


B. INVERSE DYNAMICS 





Inverse dynamics computes the states of the 
instantaneous position of the vessel along the virtual arc 


of a give trajectory. This process allows for reverse 





engineering of an executed trajectory in order to determine 
the required commands to achieve this maneuver while 


maintaining the necessary boundaries of constraints. 





In order to achieve the fore mentioned process, it is 








of note that the parameters of the reference trajectory of a 
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numerical solution are calculated in WN points evenly 





distributed throughout the virtual arc, such that 
Ar=t,(N-1)" (4-1) 


Time interval between two points is calculated as in 


equation (4-6) 








we —%*,4) +(9;- 94) 


At., =t,-t,,= 
j-l j j-l 
V,+V,, 


(452) 





In order to transition from a measurement of position along 





the virtual arc to that of velocity, the kinematics for a 


navigational solution must be achieved and are shown below 








as an example. 


Given: 
x,(t), x, (t) and therefore x, (t), x,(t) and x, (t), x, (t) 


Kinematic equations: 


i, =. =Veeos ¥ (4-3) 


i, = 2 =Vsiny (4-4) 
t 


Vax +55 (4-5) 


= arctan 2 (4-6) 
x 


Convert those to the virtual domain (wher th referenc 














d 
trajectories are developed) using the virtual speed A= 
t 
, dx, dx, dt 1 
x =—+ =—1— =—Vcos¥ (4-7) 
dt dtdt A 
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gp Oe Oe Ne aay (4-8) 
dt dt dt 2 


We will set a separate 5*-order polynomial for A(r) similar 


to that of equation (3-29) 





1 O 
0 1 
0 O 
1 T; 
0 1 
0 O 

















0 0 0 0 
0 0 0 0 lay} | 
1 0 0 0 lla] |4o 
1 2 1 3 1 4 1 5 as Ay 
ea 7 eee. a = 4-9 
’ 6 2? 0 f]lat] | A, a 
1 1 A y 
T, . Ee 7 t ; t i Ay 
(a Tag} [ar] 
1 Ty Ty Ty | 





where AJ and 4; being varied parameters and 4, 4, A, and 


A, defined as 


Ay=Vor 2,=0, Ap,=V, and A, =0 (4-10) 
Next, to address the constraints imposed on the control 
parameters (or in other words to account for the controllers 
dynamics) we also need to evaluate the following derivatives 
7 xx" tx x” 
V=VAR=A x? 4x7 + 2? (4-11) 
12 12 
EX 
: xix" — x"x! 
P= P71 = 1722 cos’ ¥ (4-12) 
1 
C. OPTIMIZATION 





In order to capitalize on the ability to optimize this 


problem a general block scheme approach is taken in order to 


accomplish this goal. 
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Define boundary conditions 
Choose reference functions 


Compute the coefficients of reference functions 


Guess on initial values of variable parameters 





Integrate speed and find states and controls 
along the virtual arc via inverse dynamics 





+ icrance 


Change the values of variable parameters 


Figure 10. General Block-Scheme 









In short, Figure 10 helps to explain that this is a 
cycle and the values of each parameter are weighted to 


ensur th rrors from the minimization function are 





continuously calculated until they are approximately zero, 





which will allow the UNREP rendezvous to be accomplished. 





Perforamance Index = » Cost Functions + » Weight x Penalty (4-13) 
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Therefore the computational procedure, shown on Figure 10 


looks as follows. We start from guessing on varied 
. " " ” " ”" " 
parameters, which are Xp, Xr Myer Apr Apr Ay and t,. Then 


using the given initial and final conditions, which are 


shown in Table 4 below, 











Initial Initial Final Final 


Position x, | Position x2 |Position x; /Position x) 





0 0 5000 yds 5000 yds 




















Table 3. Initial & Final Positions 








one can compute the coefficients of the reference functions, 


which are based off of algebraic polynomials of degree “n” 





with the virtual arc “rt” as the argument. This allows for 


independent optimization of the velocity history along the 





trajectory. Next the coefficients of th referenc 


functions are determined using MATLAB. The user then makes 








an initial guess on the values of the variable parameters. 


Using inverse dynamics, the states and controls along 


the virtual arc are determined. From this point the cost 








function and penalties are computed and if these items are 








found to within tolerance then the computation stops. If 





the tolerances are not met, the values of the variable 
parameters are changed and the states and controls along the 


virtual arc are recalculated and the cycle repeats as shown 





in Figure 10. 





Finally we estimate the performance index and evaluate 


penalties, the penalty weights are shown in Table 4. 
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Penalty Parameter Weights 














Time 10-3 
Sway Velocity 10 
Speed 10 











Table 4. Penalty Weights 


Additionally, the penalties are calculated by the equations 





shown below. 


Time Penalty = (total time - Ty (4-14) 





Speed Penalty = max([0,speed - speed,,,. ])” (4-15) 
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V. RESULTS AND CONCLUSIONS 


A. RESULTS 


Using the technique described throughout this thesis, 
the following results were achieved for three Specialized 
Command and Control Crafts (SCCC) conducting an Underway 


Replenishment at Sea with a standard supply vessel. 





For the three SCCCs, the initial and final velocities 
are the same, additionally for two out of three of the 


SCCCs, the initial angles ar th same, however all three 











have the same final angles. The change in one of the 


initial angles was done to allow for a better visual image 





of the scenario. For these results the times of arrival 





were all the same, however to accommodate the supply vessel 


all SCCCs would have varying times of arrival in order for 





the supply ship to connect one vessel then move to the next 








one and so on. Additionally, once along side, the ability 





for the real time updating at a rate of (1Hz) would allow 








the supply vessel to alter course while all SCCCS maintain 





their stations. 
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Figure 12. 
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Figure 13. 
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B. CONCLUSIONS 


This model using the fore mentioned design method 
allows the problem boundary conditions to be satisfied 
beforehand; eliminates wild trajectories which decreases 
required computer computing times; allows for real time 
updating (1Hz) of the required outputs due to the speed at 
which calculations can be done; and allows for the 
implementation of multiple agents for collision—-free 


motions. 


Additionally benefits of this technology would be that 





it allows for the incorporation of a heads up display which 
could use standard maneuvers to build a database of near- 


optimal trajectories calculated beforehand. These near- 
45 


optimal trajectories could all for the officer of the deck 








on board Naval vessels to not just mentally visualize the 





required commands, but also allow them to actually see a 


Simulation of where them will end up, thus adding to the 





overall situational awareness. In principle, if this 











onboard computer was capable of doing the required updates 


often enough, the need for a traditional feedback controller 





would be unnecessary. The computer would then be capable of 





continuously regenerating from the vessels actual position 





instead of fighting the disturbance of the trajectories. 
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APPENDIX. MATLAB CODE 


A. STARTMESCCC .M 


Q 


& This is a main script 
clear all, close all, clc 
syms tf x0 xpO xppO xf xpf xppf 





6% Setting the boudary conditions 
global posxi vxi posxf vxf 

global posyi vyi posyf vyf 

vi = 5*0.5144; vf = 13*0.5144; 
































posxi = —50; posyi = 60; % initial position, m 
posxf = 4650; posyf = 2250; % final position, m 
anglei = 0; anglef=0; 

vx = vi*cos(anglei); vyi = vi*sin(anglei); % components of initial 
velocity 

vxt = vf*cos(anglef); vyf = vf*sin(anglef); % components of final 
velocity 





6% Defining optimization problem 

















global Cost_T Fine_Speed Fine_Accel Fine_YawRate Psi_dot_max v_max a_max 
T 

global wv wvd wPsid 

R=14; % lateral seperation distance, m 





v_max=40*0.5144; 
a_max=10.7; 


ole 


maximum speed, m/s/s 
maximum acceleration, m/s/s 


ol? 




















6T = 2*sqrt ((posxf-posxi) *2+(posyf—-posyi)*2)/(vitvf)/2; % 
predetermined time of arrival, s 

T = 550; 

Psi_dot_max = 10*pi/180; % maximum yaw rate, rad/s 

wy = 10; % weighting coefficient for speed 

wvd = 10; % weighting coefficient for accelaration 


wPsid = 10; 


ol? 


weighting coefficient for yaw rate 

6% Guessing on the varied parameters 

guess (1)=0.02*sqrt ( (posxf-posxi) *2+ (posyf-posyi)*2); % virtual arc 
length 
































guess (2)=rand*0.01; % initial acceleration in x 
guess (3) =rand*-0.0001; % initial acceleration in y 
guess (4) =rand*-0.0001; % final acceleration in x 

guess (5)=rand*0.001; & final acceleration in y 

guess (6) =rand*-0.0001; § initial acceleration in lambda 
guess (7)=rand*0.0001; % final acceleration in lambda 


6% Defining coefficients (units are commented) 


sDefineSCCC 

6% Compute Coefficients 

global aM 

A=[1 0 0 0 0 0; 
O01 0 0 0 0; 
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00 1 0 0 Os 
1 tf t£%2/2 t£%3/6 tf*4/12 t£%5/20; 
Of. “tt t£*2/2 t£%3/3 t£*4/4; 
00 1 tee bE? ESS]? 
b=[x0 xpO xppO xf xpf xppf]'; 
a=A\b; 


a=collect(a,tf); 
M=length (a); 





6% Calling the optimization routine 

opt=optimset ('Display', 'iter', 'TolX',1le-4, 'TolFun',1le-4, 'MaxIter',10); 
[guess_opt, fval,exitflag]=fminsearch('Trajectory',guess,opt) ; 
sTrajectory (guess) ; 


6% Displaying cost function and penalties 





fprintf('Cost function 2 ©64S6.20\n', Cost_T) 
fprintf£(' Penalty in speed : %6.2g\n',Fine_Speed) 
fprintf£(' Penalty in acceleration : %6.2g\n',Fine_Accel) 
fprintf(' Penalty in yaw rate : %6.2g\n\n',Fine_YawRate) 


6% Displaying optimal parameters 
Sguess_opt=guess; 
fprintf('Arc lenght = %6.2f\n',guess_opt(1)) 








fprintf(' initial aces! final accel\n') 
fprintf('along x : $6.2e $6.2e\n', [guess_opt (2) 
guess_opt (3) ]) 

fprintf('along y : $6.2e $6.2e\n', [guess_opt (4) 
guess_opt (5) ]) 

fprintf('in lambda: $6.2¢e %$6.2e\n', [guess_opt (6) 


guess_opt (7) ]) 


6% Plotting the results 
PlotResults 


B. TRAJECTORY .M 


function PlI=trajectory (guess) 

6% This function computes states and controls for the current guess 
lobal aM 

obal tf x0 xp0O xppO xf xpf xppf 

yms tf x0 xp0 xpp0O xf xpf xppf 

obal posxi vxi posxf vxf 

obal posyi vyi posyf vyf 

obal t xl x2 v Psi tau lam v_dot Psi_dot 




















QQqQgaugqea 





6% Current values of varied parameters 









































tauf =guess(1); % virtual arc length 

accxi=guess (2); % initial acceleration in x 
accyi=guess (3); % initial acceleration in y 
accxf=guess (4); & final acceleration in x 
accyf=guess (5); % final acceleration in y 
accli=guess (6); 6 initial acceleration in lambda 
acclf=guess (7); & final acceleration in lambda 
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6% Defining coordinates in N nodes in the virtual domain 
al = subs (a, {('sx0", "spo? "xpp0!, 'xz", "xpk", "epee! ,. EE" by eas 


{posxi,vxi,accxi,posxf,vxf,accxf,tauf}); 
az = subs (a, { "20", 'xpa", *spe0", "xi", spe’, Vxppet, "ER" },e44 


{posyi,vyi,accyi,posyf,vyf,accyf,tauf}); 
a3: = subs (a; {('x0", 'xpo?; 'xpp0", 'xrt, xp", eppe!, TEE! } yaa 








{1,0,accli,1,0,acclf,tauft}); 

axl=diag([1,1,1/2,1/6,1/12,1/20]) *al; 
ax2=diag([1,1,1/2,1/6,1/12,1/20]) *a2; 
ax3=diag([1,1,1/2,1/6,1/12,1/20]) *a3; 





tau=linspace(0,tauf) ; 




















6% Defining coordinates' derivatives in N nodes in the virtual domain 





cxl_prime = cx1l.*[5:-1:0]*eye(6,5); 
cx2_prime = cx2.*[5:-1:0]*eye(6,5); 
cx3_prime = cx3.*[5:-1:0]*eye(6,5); 
xl_prime = polyval(cxl_prime,tau) ; 


x2_prime = polyval 
lam_prime = polyval 


cx2_prime,tau) ; 
cx3_prime,tau); 














( 
( 


6% Defining coordinates' second-order derivatives in N nodes 
exl_2prime=cxl_prime.*[4:-1:0]*eye(5,4); 
ex2_2prime=cxl_prime.*[4:-1:0]*eye(5,4); 
x1_2prime=polyval (exl_2prime, tau) ; 
x2_2prime=polyval (ex2_2prime,tau) ; 























N=length (x1); 





del_tau tauf/ (N-1); 
t (1) = 0; % time 
v(1) = sqrt (vxi*2+vyi*2); % initial speed, m/s 
Psi = atan2(x2_prime, xl_prime) ; % heading, rad 
for j=2:N 
sq = sqrt ((xl_prime(j))%*2+(x2_prime(j))%2); 
v(j) = lam(j)*sq; % speed, m/s 
dt = 2*sqrt ( (x1 (3) -x1 (5-1) ) *2+ (x2 (3) -x2 (3-1) ) *%2)/ (v (4) +v (5-1) ) 5 
t(j) = t(j-1)+dt; 
v_dot (4) = lam_prime(j)*sqt... 
lam(4j) *2* ((xl_prime (j) *x1_2prime (j))+(x2_prime (j) *xl_2prime(j)))/sq; 
Psi(j) = atan2(x2_prime(j),xl_prime(j)); 
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Psi_dot(j) = (xl_prime(j)*x2_2prime(j)-x1l_2prime(j)*x2_prime(j))... 


/x1l_prime(j)*2*cos (Psi(4j))%2; 


end 
PI = PerformancelIndex; 
reluEn 


function PI=PerformanceIndex 
% This function computes the combined performance index 


% 


P 


Q@Haa 


lobal 


t xl x2 v Psi tau lam v_dot Psi_dot 





lobal 











obal 





Cost_T Fine_Speed Fine_Accel Fine _YawRate Psi_dot_max v_max a_max 


wv wvd wPsid 





Cost_T = (t(end)-T)%2; 
Fine Speed = max([0,abs (v) -v_max])*2; 


Fine_Accel 
Fine_YawRate 


max([0, (abs (v_dot) -a_max) ])%*2; 
max([0, (abs (Psi_dot) —Psi_dot_max) ]) *2; 


( 
Fine_A = max(0,R-min(sqrt ((xl-posxfr) .*%2+(x2- 
posyEr) <2) ))*23 
I = Cost_T + wv*Fine_Speed + wvd*Fine_Accel + wPsid*Fine_YawRate; 
return 


Cc. 


PLOTRESULTS .M 


is script plots the results of optimization 


t xl x2 v Psi tau lam v_dot Psi_dot 














Cost_T Fine_Speed Fine_Accel Fine _YawRate Psi_dot_max v_max a_max 





%% Bird-eye view 
figure ('Name','2D View') 
title('Optimal Trajectory') 


ro) 
oO 





lot (x1,x2) 
label('x_1, m'), ylabel('x_2, m') 
xis equal, grid on, hold on 

lot (x1 (1),x2(1),'r0O") 

lot (x1 (end) ,x2 (end), 'rO') 











% Plotting time histories 





figure ('Name', 'States') 
ubplot (2,1,1) 


s 
Pp 


x] 


Pp 





abel 


lot(t,v), hold 


('Time, s'), ylabel('Speed, m/s') 





lot 


— 


0 t(end)],v_max*[1 1],'r--') 


subplot (2,1,2) 


Pp 


x] 





abel 


lot (t,Psi*180/pi) 


A 





('Time, s'), ylabel('Heading, *o') 





figure ('Name', 'Controls') 
subplot (2,1,1) 


Pp 
x 
p 
p 





lot (t,v_dot), hold 

label ('Time, s'), ylabel('Acceleration, m/s*2"') 
lot ({[0 t(end)],a_max*[1 1], 'r--') 

lot ([0 t(end)],-a_max*[1 1],'r--') 


subplot (2,1,2) 
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plot (t,Psi_dot*180/pi), hold 

xlabel('Time, s'), ylabel('Yaw rate, “o/s') 

plot ([0 t (end) ],Psi_dot_max*[1 1]*180/pi,'r--') 
plot ([0 t (end) ],-Psi_dot_max*[1 1]*180/pi,'r--') 





figure('Name', 'Virtual Domain Parameters') 

subplot (2,1,1) 

lot (t, tau) 

label ('Time, s'), ylabel('\tau') 

ubplot (2,1,2) 
) 
m 











lot (t, lam 
label ('Ti 








xX OH KO 























e, 8'), ylabel ("\lambda, s*{-1}") 
D. TEST .M 


close all 

SupplyShip (4700,2200,0.2) 
SCCCship (50,0,1) 
SCCCship (40, 90,1) 
SCCCship (-50, 60,1) 


E. SUPPLYSHIP .M 


function SupplyShip (X,Y,Psi) 
% X,Y - the coordinates of the ship's center in the local tangent plane 
in meters 
% Psi - the orientation of the ship in the local tangent plane in 
radians 
SCCCx=[0 250 700 754.6 700 250 0]*0.3048; 
SCCCy=[100 107 107 53.7 0 0 7)*0.3048; 
SCCCx=SCCCx-mean (SCCCx) ; 
SCCCy=SCCCy-mean (SCCCy) ; 
len=sqrt (SCCCx.*2+SCCCy.%*2); 
ang=atan2 (SCCCy, SCCCx) ; 
SCCCx=Xt+len.*cos(Psi-ang) ; 
SCCCy=Ytlen.*sin(Psi-ang) ; 
patch (SCCCx, SCCCy, 'c') 
axis equal 











F. SCCCSHIP .M 


function SCCCship(X,Y,Psi) 
% X,Y - the coordinates of the ship's center in the local tangent plane 
in meters 
% Psi - the orientation of the ship in the local tangent plane in 
radians 
SCCCx=[0 130 135 130 0]*0.3048; 
SCCCy=[68 68 34 0 0]*0.3048; 
SCCCx=SCCCx-mean (SCCCx) ; 
SCCCy=SCCCy-mean (SCCCy) ; 
len=sqrt (SCCCx.*2+SCCCy.%*2); 
ang=atan2 (SCCCy, SCCCx) ; 
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SCCCX=X+4 





SCCCy=Y4 


len.*cos (Psi-ang) ; 





len.*sin(Psi-ang) ; 


patch (SCCCx, SCCCy, 'm') 


axis equal 
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